library(dplyr)
library(tidyr)
library(ggplot2)

library(ggrepel)
library(ProbBayes)

curve(dnorm(x, mean = 0, sd = 1), 
      from = -10, to = 10, 
      col = "red", lwd = 2, label = "Normal Distribution",)

curve(dcauchy(x, location = 0, scale = 1), 
      from = -10, to = 10, 
      main = "Cauchy Distribution", 
      xlab = "x", ylab = "Density", 
      col = "blue", lwd = 2, add= TRUE, label = "Cauchy Distribution")

# legend("topright", legend = c("Normal Distribution", "Cauchy Distribution"), 
#        col = c("red", "blue"))


p_location <- c(0,0,1,0,0,0)

p_prob <- matrix(c(.5, .5, 0, 0, 0, 0,
                   .25, .5, .25, 0, 0, 0,
                   0, .25, .5, .25, 0, 0,
                   0, 0, .25, .5, .25, 0,
                   0, 0, 0, .25, .5, .25,
                   0, 0, 0, 0, .5,.5), nrow = 6, ncol = 6, byrow=TRUE)

print(p_location %*% p_prob, digit =5)

print(p_location %*% p_prob %*% p_prob %*% p_prob %*% p_prob %*% p_prob, digit = 5)


for (xi in 1:2000){
  p_location <- p_location %*% p_prob
}

print(p_location, digit = 5)

pm <- diag(rep(1,6))

for(xi in 1:2000){
  pm <- pm %*% p_prob
}

print(pm, digit = 5)


s <- vector("numeric", 10000)
s[1] <- 4

for (i in 2:10000) {
  s[i] <- sample(1:6, size=1, prob = p_prob[s[i-1],])
  
}


s_data <- data.frame(Iteration = 1:10000, 
                Location = s)

s_data %>% mutate(
  L1 = (Location == 1),
  L2 = (Location == 2),
  L3 = (Location == 3),
  L4 = (Location == 4),
  L5 = (Location == 5),
  L6 = (Location == 6)
) %>%
  mutate(Proportion1 = cumsum(L1) / Iteration,
         Proportion2 = cumsum(L2) / Iteration,
         Proportion3 = cumsum(L3) / Iteration,
         Proportion4 = cumsum(L4) / Iteration,
         Proportion5 = cumsum(L5) / Iteration,
         Proportion6 = cumsum(L6) / Iteration) %>%
  select(Iteration, 
         Proportion1, Proportion2, Proportion3, 
         Proportion4, Proportion5, Proportion6) -> s_data_1

gather(s_data_1, Outcome, Probability, -Iteration) -> S2

ggplot(S2, aes(x = Iteration, y = Probability, color = Outcome)) +
  geom_line() +
  facet_wrap( ~ Outcome, ncol=3) + 
  ylim(0, .4) +
  ylab("Relative Frequency") +
  theme(text= element_text(size = 18)) +
  scale_color_manual(values = c("red", "blue", "green", "purple", "orange", "brown")) +
  scale_x_continuous(breaks = seq(0, 10000, by = 1000)) +
  theme_minimal() 
  

pd <- function(x){
  values <- c(5,10,4,4,20,20,12,5)
  ifelse (x %in% 1:length(values), values[x], 0)
}

pd(1:8)


prob_dist <- data.frame(x=1:8,
                        prob=pd(1:8))

prob_plot(prob_dist, Color= 'blue') +
  theme(text=element_text(size=18)) +
  ylab("Probability")

random_walk <- function(pd, start, num_steps){
  y <- rep(0, num_steps)
  
  current <- start
  
  for (xi in 1:num_steps) {
    candidate <- current + sample(c(-1, 1), size = 1)
    prob <- pd(candidate) / pd(current)
    if (runif(1) < prob) {
      current <- candidate
    }
    y[xi] <- current
  }
  
  return (y)
}

replicate(10, runif(1))


out <- random_walk(pd, start = 4, num_steps = 10000)

data.frame(out) %>%
  group_by(out) %>%
  summarise(N=n(), Prob=N/10000) -> s

prob_dist %>% mutate(Prob = prob /sum(prob)) -> prob_dist


df <- rbind(data.frame(x=prob_dist$x,
                       Prob=prob_dist$Prob,
                       Type='actual'),
            data.frame(x=s$out, Prob=s$Prob, Type='Simulated')
)

ggplot(df, aes(x=x, y=Prob, fill=Type)) +
  geom_bar(stat='identity', position=position_dodge(), color='black') +
  theme(text=element_text(size=18)) +
  ylab("Probability") +
  xlab("x") +
  ggtitle("Actual vs Simulated Probability Distribution") +
  theme_minimal()
